    % section 6.4.1 'robustness' / EZ 
        % this file shows the differences in CE for GHHW and GLTHI 
        % using Epstein Zin preferences
    
    clear;
    load('./results/allstates_main.mat');
    load('./results/model_parameters_main.mat');
    load('./results/capHtmod_c.mat');
    
    addpath functions;

    rho = model_parameters.rho;
    MU = model_parameters.MU;
    CAPH = model_parameters.CAPH;
    gamma = model_parameters.gamma;
    iX = simulation_output.iX;
    Nstates = 7;
    Ppaid_N = simulation_output.Ppaid_N;
    Wio_abi = simulation_output.Wio_abi;
    Wio_rea = simulation_output.Wio_rea;

    maxreps = 100;
    tolFP=0.005;
    T = 70;
    
    
   r = 1/rho-1;
    
   rng(1234567,'twister')

   prob_start_fixed_abi = capHtmod_c;
   psi_grid = [1.1:0.1:4.1];
   gamma1 = 0.000075;
   gamma2 = 0.0004;
   gamma3 = 0.0008;
  
   CE_PHI_CARA_abi_all1 = arrayfun(@(psi)EZ_welfare_PHI_v4(Wio_abi*1000,CAPH,MU*1000,maxreps,tolFP,rho,psi,gamma1,@invG_cara,@G_cara,prob_start_fixed_abi),psi_grid);
   CE_HHW_CARA_abi_all1 = arrayfun(@(psi)EZ_welfare_HHW_v4(Wio_abi*1000,CAPH,MU*1000,maxreps,tolFP,rho,psi,gamma1,@invG_cara,@G_cara,prob_start_fixed_abi),psi_grid);
   DELTACE_EZ1 =  (CE_HHW_CARA_abi_all1-CE_PHI_CARA_abi_all1)./CE_HHW_CARA_abi_all1;

   CE_PHI_CARA_abi_all2 = arrayfun(@(psi)EZ_welfare_PHI_v4(Wio_abi*1000,CAPH,MU*1000,maxreps,tolFP,rho,psi,gamma2,@invG_cara,@G_cara,prob_start_fixed_abi),psi_grid);
   CE_HHW_CARA_abi_all2 = arrayfun(@(psi)EZ_welfare_HHW_v4(Wio_abi*1000,CAPH,MU*1000,maxreps,tolFP,rho,psi,gamma2,@invG_cara,@G_cara,prob_start_fixed_abi),psi_grid);
   DELTACE_EZ2 =  (CE_HHW_CARA_abi_all2-CE_PHI_CARA_abi_all2)./CE_HHW_CARA_abi_all2;
   
   CE_PHI_CARA_abi_all3 = arrayfun(@(psi)EZ_welfare_PHI_v4(Wio_abi*1000,CAPH,MU*1000,maxreps,tolFP,rho,psi,gamma3,@invG_cara,@G_cara,prob_start_fixed_abi),psi_grid);
   CE_HHW_CARA_abi_all3 = arrayfun(@(psi)EZ_welfare_HHW_v4(Wio_abi*1000,CAPH,MU*1000,maxreps,tolFP,rho,psi,gamma3,@invG_cara,@G_cara,prob_start_fixed_abi),psi_grid);
   DELTACE_EZ3 =  (CE_HHW_CARA_abi_all3-CE_PHI_CARA_abi_all3)./CE_HHW_CARA_abi_all3;

   F_ez13 = [psi_grid',DELTACE_EZ1',DELTACE_EZ2',DELTACE_EZ3'];
   writetable(array2table(F_ez13),'./results/F_ez13.xlsx')    
   
   % Ed 10
    prob_start_fixed_rea = capHtmod_c;


   CE_PHI_CARA_rea_all1 = arrayfun(@(psi)EZ_welfare_PHI_v4(Wio_rea*1000,CAPH,MU*1000,maxreps,tolFP,rho,psi,gamma1,@invG_cara,@G_cara,prob_start_fixed_rea),psi_grid);
   CE_HHW_CARA_rea_all1 = arrayfun(@(psi)EZ_welfare_HHW_v4(Wio_rea*1000,CAPH,MU*1000,maxreps,tolFP,rho,psi,gamma1,@invG_cara,@G_cara,prob_start_fixed_rea),psi_grid);
   DELTACE_EZ1 =  (CE_HHW_CARA_rea_all1-CE_PHI_CARA_rea_all1)./CE_HHW_CARA_rea_all1;

   CE_PHI_CARA_rea_all2 = arrayfun(@(psi)EZ_welfare_PHI_v4(Wio_rea*1000,CAPH,MU*1000,maxreps,tolFP,rho,psi,gamma2,@invG_cara,@G_cara,prob_start_fixed_rea),psi_grid);
   CE_HHW_CARA_rea_all2 = arrayfun(@(psi)EZ_welfare_HHW_v4(Wio_rea*1000,CAPH,MU*1000,maxreps,tolFP,rho,psi,gamma2,@invG_cara,@G_cara,prob_start_fixed_rea),psi_grid);
   DELTACE_EZ2 =  (CE_HHW_CARA_rea_all2-CE_PHI_CARA_rea_all2)./CE_HHW_CARA_rea_all2;
   
   CE_PHI_CARA_rea_all3 = arrayfun(@(psi)EZ_welfare_PHI_v4(Wio_rea*1000,CAPH,MU*1000,maxreps,tolFP,rho,psi,gamma3,@invG_cara,@G_cara,prob_start_fixed_rea),psi_grid);
   CE_HHW_CARA_rea_all3 = arrayfun(@(psi)EZ_welfare_HHW_v4(Wio_rea*1000,CAPH,MU*1000,maxreps,tolFP,rho,psi,gamma3,@invG_cara,@G_cara,prob_start_fixed_rea),psi_grid);
   DELTACE_EZ3 =  (CE_HHW_CARA_rea_all3-CE_PHI_CARA_rea_all3)./CE_HHW_CARA_rea_all3;

  
   F_ez10 = [psi_grid',DELTACE_EZ1',DELTACE_EZ2',DELTACE_EZ3'];
   writetable(array2table(F_ez10),'./results/F_ez10.xlsx')    


         
         

